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The recent discovery of tunable Dzyaloshinskii-Moriya interactions in layered magnetic materials 
with perpendicular magnetic anisotropy makes them promising candidates for stabilization and 
manipulation of skyrmions at elevated temperatures. In this article, we use Monte Carlo simulations 
to investigate the robustness of skyrmions in these materials against thermal fluctuations and finite- 
size effects. We find that in confined geometries and at finite temperatures skyrmions are present in 
a large part of the phase diagram. Moreover, we find that the confined geometry favors the skyrmion 
over the spiral phase when compared to infinitely large systems. Upon tuning the magnetic field 
through the skyrmion phase, the system undergoes a cascade of transitions in the magnetic structure 
through states of different number of skyrmions, elongated and half skyrmions, and spiral states. We 
consider how quantum and thermal fluctuations lift the degeneracies that occur at these transitions, 
and find that states with more skyrmions are typically favored by fluctuations over states with less 
skyrmions. Finally, we comment on electrical detection of the various phases through the topological 
and anomalous Hall effects. 


I. INTRODUCTION 


A skyrmion is a certain type of topological field config¬ 
uration which was first introduced in particle physics. It 
corresponds to a classical stationary solution of the equa¬ 
tions of motion with which one can associate a topological 
invariant and describes the emergence of a discrete par¬ 
ticle from a continuous field^i More recently, skyrmions 
have been considered in quantum Hall devices^ Bose- 
Einstein condensatesand liquid crystalsMag¬ 
netic skyrmions were predicted^ and recently observed 
in bulk materials like MnSi and Cu20Se03 at low 
temperaturesIn these materials, bulk inversion sym¬ 
metry is broken which allows for nonzero Dzyaloshinskii- 
Moriya (DM) interaction o^^d^ that lead to helical, con¬ 
ical, and skyrmionic spin textures depending on the ap¬ 
plied magnetic field^^ It was proposed that skyrmions 
are promising candidates for encoding binary data that 
allow for high-density and low power consumption mag¬ 
netic memories due to low critical currents for skyrmion 
motion and their inherent topological stability^i^”— In 
addition to fundamental interest concerning the interplay 
between topology, geometry and fluctuations, it is there¬ 
fore relevant to understand the behavior of skyrmions at 
high temperature and in confined thin-film geometries for 
realizing spintronic devices. 

Motivated by recent experiments on domain wall mo¬ 
tion in magnetic thin films with perpendicular magnetic 
anisotropy (PMA materials) that point to sizeable and 
tunable DM interactionsi^— and anisotropy, we use 
Monte Carlo simulations to investigate the robustness of 
skyrmions in these systems in confined wire-like geome¬ 
tries against thermal fluctuations. Our main results are: 
i) The phase diagram in Fig.[T]that shows which spin tex¬ 
tures are to be expected at a certain anisotropy strength 


and applied magnetic field at nonzero temperature. We 
find that the confining geometry extends (with respect 
to systems in the thermodynamic limit) the skyrmion 
phase at the expense of the spiral phase, ii) The cascade 
of transitions between different magnetic structures that 
the system undergoes upon lowering the magnetic field 
through the skyrmion phase, a typical example of which 
is shown in Fig. [3l iii) The temperature dependence of 
relative probabilities for occurence of different skyrmions 
configurations that are degenerate at zero temperature, 
shown in Fig. [6l We find that at moderate tempera¬ 
ture configurations with more skyrmions are typically en- 
tropically favored over configuration with less skyrmions. 
Moreover, we also consider quantum fluctuations and find 
that these also favor configurations with higher skyrmion 
number. 

Skyrmions are predicted to occur in several varieties^ 
two of which have by now been experimentally identified. 
One of these is the Neel (somtimes also called “hedgehog”) 
skyrmion, i.e., a skyrmion in which the magnetization 
points radially outward from the skyrmion center. The 
other type of skyrmion (where the magnetization is per¬ 
pendicular to radii pointing outward from the skyrmion 
center) is called a Bloch skyrmion. (See Fig. I of Ref. 
for an illustration.) Which of the two types is favored 
depends on which type of DM interactions are present as 
dictated by the crystal and/or structural symmetry. In 
the long-wavelength limit, the DM interactions yield a 
contribution to the energy density that is a certain com¬ 
bination of so-called Lifshitz invariants, i.e., antisymmet¬ 
ric terms of the form 


dSj ^ 

* Qj, j Qj, 


with Si the i-th cartesian component of the spin and 
r = x,y oi z, di spatial direction. In MnSi, one of the 
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most-studied skyrmion materials, the DM interactions 
give a contribution proportional to 

S-(VxS), (2) 

which is straightforwardly shown to be a combination 
of terms of the form as in Eq. o, and favors Bloch 
skyrmions. For the PM A materials that are the focus 
of this work, the DM interactions are interface-induced 
(see also Ref. [23 ) and stabilize Neel skyrmions. They are 
proportional to the expression 

(z.S)(V.S)-(S.V)(z.S) , (3) 

and in this particular form are shown to explicitly de¬ 
pend on the symmetry-breaking direction z which de¬ 
notes the normal to the interface. We note, however, that 
the above form of DM interactions also arises in crystals 
of symmetry class Cnv and therefore that Neel skyrmions 
are stabilized by bulk DM interactions in some materials. 
This was recently observed in the magnetic semiconduc¬ 
tor GaV 4 S 8 '^ In this material there are no Lifshitz in¬ 
variants in the 2 ;-direction and a conical phase is therefore 
not present. Our results therefore also apply to this case, 
and, particularly when grown in thin-film form on differ¬ 
ent substrates, the magnetic anisotropy of this material 
may be tuned between easy axis and easy plane which 
enables experimentally exploring the full phase diagram 
in Fig. [TJ 

Previous theoretical studies have focused on the 
skyrmion phase diagram of infinite systems at zero 
temperature^i^ and nonzero temperature in two and 
three dimensionsThe first experimental results on 
bulk materials^— have been extended to confined ge¬ 
ometries, such as thin films, e.g. of MnSi^ and FeGe4^ 
Transitions between states with a different number of 
skyrmions and other (non-skyrmion) magnetic configu¬ 
rations as a function of field have been discussed experi¬ 
mentally and theoretically (at zero temperature) for thin 
films of MnSi (and thus for bulk DM interactions leading 
to Bloch skyrmions) in Refs. lail and [32I. Evidence for 
such a cascade of transitions in magnetoresistance mea¬ 
surements was very recently discussed in Ref. [ssj. In these 
works, the Bloch skyrmion core lines (and field) lie in 
the plane of the thin film or along the wire direction due 
to strong easy-plane anisotropy induced by tensile strain 
from the substrate. As a result, only in-plane fields sta¬ 
bilize skyrmions whereas fields perpendicular to the thin 
film (and thus anisotropy plane) lead to conical phases. 
In the two-dimensional configuration that we consider, 
however, (see Fig. [3|) the external field is perpendicu¬ 
lar to the thin film and skyrmions may be observed for 
both easy-plane and easy-axis anistropy. For the form of 
the DM interactions considered here [Eq. ([3])] the conical 
phase is absent and the easy-plane anisotropy leads to a 
large external-field range over which skyrmions are stable 
(see Fig. [T]). Moreover, the confined geometry stabilizes 
skyrmions over spirals as the skyrmion phase becomes 
larger in the wire geometry with respect to the case of 
an infinite system. 


Very recently, Du et al reported on real-space observa¬ 
tion of a cascade of transitions in the magnetic structure 
of a FeGe wire^^ Although this system stabilizes Bloch 
skyrmions, this work gives experimental corroboration of 
many of our findings, such as the existence of elongated 
skyrmions, the creation of skyrmions from edge states, 
and the enhanced stability of skyrmions in the confining 
geometry. Moreover, these authors also observe that spi¬ 
rals orient themselves with their wave vector along the 
edge of the wires (see also Ref. [35|) , which is in line with 
our results as well. 

Finally, we note for completeness that single Bloch 
skyrmions in nanodisks were considered in Ref. [sH with¬ 
out taking into account anisotropy and that the effects 
of thermal fluctuations and confining geometry on single- 

e rmion dynamics was considered in Refs. [^, [^, and 

The remainder of this article is organized as follows. 
In Sec. mi we discuss our model hamiltonian and the al¬ 
gorithm used in our simulations. Using these simulations 
we then construct the anisotropy-external field phase dia¬ 
gram for wire-like systems at moderate temperature and 
discuss the various magnetic phases and cascade of tran¬ 
sitions between them in Sec. mil Hereafter we focus on 
points in the phase diagram where degeneracies occur 
at zero temperature and investigate how fluctuations lift 
these degeneracies. We end with a conclusion and out¬ 
look, and briefly discuss electrical detection of the various 
magnetic configurations. 



Figure 1. The anisotropy-external field phase diagram for 
a wire geometry with L = 16 and p = 8 at ksT/J = 0.5 
(solid line) and zero temperature (dashed line). The dot¬ 
ted line corresponds to the case of an infinite system at zero 
temperature.— These lines encircle the region where the wind¬ 
ing number is larger than one half. The susceptiblity of the 
winding number at finite temperature Xw is also indicated. 
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II. MODEL AND SIMULATIONS 

We consider Heisenberg spins Sr of unit length, where 
r denotes the location on a two-dimensional square lattice 
in the x-^-plane. The spins interact via a ferromagnetic 
Heisenberg coupling J > 0, Dzyaloshinskii-Moriya cou¬ 
pling D, and are subject to an anisotropy term K and 
an external magnetic field B. The effective hamiltonian 
is given by 

H = — J Sr • (Sr+x + ^r+y) 

r 

+ K^(Sr-zf -B-^Sr (4) 

r r 

-D y] (Sr X Sr+x • y - Sr X Sr+y ' x) . 

r 

The above form of the DM interactions (i.e., the term in 
the hamiltonian proportional to D) corresponds to the 
discretized version of Eq. (|3|) and is therefore to lowest 
order in nearest-neighbor coupling appropriate for the 
PMA materials of interest to us here, i.e., for the sit¬ 
uation that the DM interactions arise due to inversion 
asymmetry induced by the presence of an interface. In 
the hamiltonian we have neglected dipole-dipole inter¬ 
actions. This is appropriate in the limit of strong DM 
interactions that lead to small skyrmions and where the 
dipolar field will only renormalize parameters such as the 
anisotropy. We note that throughout this article we con¬ 
sider a wire geometry, i.e., a two-dimensional system with 
periodic boundary conditions in one direction and open 
boundary conditions in the other direction. 

For simplicity we assume B = Hz perpendicular to the 
x-^-plane and take H > 0 without loss of generality. The 
dimensionless parameters ksT/J = l/(ydJ), with ksT 
the thermal energy, and KJjD^ determine the 

state of the system. An important length scale is the 
pitch length p (in units of the lattice constant) that de¬ 
termines the periodicity of magnetic textures that arise 
due to competition between DM interaction D and ex¬ 
change J. In case of a spiral state, for example, the 
pitch length will be the period of the spirals. In the 
absence of anisotropy, A = 0, the pitch length is given 
by2i D/J = tan(27r/p) and we will use this definition 
throughout. This relation can be used to coarse-grain 
the system by keeping the ratio between the dimensions 
of the system and the pitch length constant while chang¬ 
ing discretization. 

We use classical Monte Carlo simulations to sample 
phase space at finite temperatures T. A typical simula¬ 
tion starts with a completely randomized square lattice 
of spins with Lx L sites (with one open and one periodic 
boundary) at fixed parameters BJ/D^ and KJjD^ at a 
scaled temperature ksT/ J = 10.0 far above the critical 
temperature of the system. After several lattice updates 
the temperature is lowered to a fraction of 0.95 of the 
last temperature until a temperature ksT/J = 0.01 is 
reached. At each temperature data is collected and 100 


simulations are independently run per set of parameters. 

At every temperature the lattice is updated using the 
Metropolis algorithm followed by a cluster-flipping algo¬ 
rithm. The Metropolis algorithm consists of choosing a 
spin at random from the lattice and proposing a new 
direction for this spin. The proposed spin direction is 
chosen uniformly from the area that is the spherical cap 
around the original direction where the maximal angle 
between the original and the proposed spin direction is 
a, and a is dynamically adjusted so that the acceptance 
probability is approximately 50 %. The cluster-flipping 
algorithm grows a cluster of spins and flips all the spins 
in the opposite direction similar to the Wolff algorithm 
for Ising spins. In both cases acceptance rates for chang¬ 
ing a spin (or cluster of spins) are based on the energy 
difference between the spin states before and after the 
move such that the system is sampled according to the 
Boltzmann distribution^^ 


BJID^ 

0.0 0.5 1.5 



(a) (b) (c) 



Figure 2. Simulations snapshots of spin systems of size L x 
L = 32 X 32 for different parameter values and temperatures 
with the parameter D/J chosen such that the pitch length 
p = 8. The periodic boundary is drawn as a red solid line and 
spin vectors are represented by colored arrows where the color 
scales linearly with the z-component of the vector to highlight 
spin textures. From top to bottom ksT/J takes on values 
5.00, 0.50, and 0.01 respectively and the system undergoes a 
transition from an unpolarized state to an ordered state in 
this temperature range. In all cases KJ/D^ is zero and from 
left to right BJ/D^ has values 0.0, 0.5, and 1.5, resulting in 
a spiral, skyrmion and polarized state, respectively, for low 
enough temperatures. 
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III. MAGNETIC PHASES AND PHASE 
DIAGRAM 

Throughout this article we use a nonstandard defini¬ 
tion of phases and phase diagram. In the systems we con¬ 
sider there can, strictly speaking, not be any thermody¬ 
namic phase transition breaking a continuous symmetry 
for two reasons: (a) we always work at finite temperature 
in a two-dimensional system where the Mermin-Wagner 
theorem forbids the breaking of a continuous symme¬ 
try even in the thermodynamic limit; (b) we explicitly 
consider finite size systems where a real thermodynamic 
phase transition is also ruled out. Nonetheless, it is pos¬ 
sible to identify phases that are distinct in their magnetic 
configurations, such as phases that may or may not have 
magnetic skyrmions. 

In order to distinguish such phases, a useful quantity 
is the winding number re, which plays the role of a topo¬ 
logical charge. It is quantized in a system with peri¬ 
odic boundary conditions and in terms of the unit vector 
n = (S)/|(S)| (we choose n instead of S to distinguish 
between the microscopic spin and the order parameter 
field) it reads 

w = ^ [ dxdyn • {dxn x dyu) , (5) 

47r J 

in the continuum limit. In the definition of n and 
throughout this article angular brackets denote the ex¬ 
pectation value in the canonical ensemble. A single 
(anti-) skyrmion contributes (minus) one to the wind¬ 
ing number. This effectively means that the number of 
skyrmions in an ordered state can be counted. From 
a set of independent simulations with identical param¬ 
eters the susceptiblity of the winding number Xw = 
((rc^) — (rc)^) J/ksT is calculated. In the latter ex¬ 
pression the expectation values are determined by us¬ 
ing Eq. m, with n replaced by S and the integral re¬ 
placed by a sum, for a given spin configuration in the 
simulations and then averaging over many such config¬ 
urations. Whenever the system undergoes a transition 
between states with different number of skyrmions this 
susceptibility will be enhanced. 

To determine the phase diagram we take the linear 
system size L equal to 16 with one periodic boundary 
condition (i.e., a wire geometry) at fixed parameter val¬ 


Finally, we note that the phase diagram at zero tem¬ 
perature was studied before for infinite system size with 
different methodsContrary to this latter work we do 
not find a tri-critical point where polarized, spiral and 
skyrmion phases meet in the easy-plane part of the phase 
diagram from our simulations. Representative spin con¬ 
figurations of the polarized state, the spiral state, and 
the skyrmion state at small temperature can be found 
in Fig. [2] (g), (h), and (i), respectively. We also note 
that for Bloch skyrmions the role of anisotropy on the 


ues BJjD^ and KJjD^. We note here that we simulate 
a square lattice for convenience throughout, and have 
not performed finite-size scaling as a function of wire 
length as we expect our findings to merely change quan¬ 
titatively for larger system size in the periodic direction. 
We slowly cooled down the system below the critical tem¬ 
perature and measured the susceptibility of the winding 
number Xw- For ksT / J = 0.5 and zero temperature we 
find the phase diagram in Fig. [H divided into an easy- 
axis {K < 0) and easy-plane {K > 0) part. Within the 
solid line in this figure, the average winding number is 
larger than one half, so that magnetic skyrmion config¬ 
urations are expected. The dashed line corresponds to 
the zero temperature case for the same wire-like confined 
geometry. The dotted line corresponds to the infinite 
system at zero temperature4^ Below the skyrmion phase 
boundary one finds a phase where spiral states rather 
than skyrmions are stabilized in the infinitely large sys¬ 
tem. At elevated temperatures {ksT/J = 0.5) and in 
the confined geometry we also find spiral states for fields 
below the skyrmion phase, albeit that more complicated 
textures also appear (depending on system size, pitch 
length, and anistropy; below we discuss the magnetic 
configurations in this region in more detail). In the con¬ 
fined geometry the skyrmion phase becomes larger (with 
respect to the infinite system) at the expense of the spiral 
phase. We attribute this to the larger ability of skyrmions 
as compared to spirals to adapt to the repulsive forces 
away from the edges of the system^^ Du et al. have 
observed this enhanced stability for Bloch skyrmions in 
FeGe nanowires 

For fields too large to stabilize skyrmions, the spins are 
uniformly polarized. In this part of the phase diagram 
and for fields larger than 2K (with K > 0) the spins are 
pointing along the field, i.e., along the ^-direction. For 
smaller fields and K > 0 the easy-plane anisotropy tilts 
the spins away from the field direction. 

The colors in the phase diagram indicate the suscep¬ 
tibility of the winding number. Depending on system 
size relative to pitch length, the magnetic configuration 
within the skyrmion phase may undergo transitions be¬ 
tween phases with different number of skyrmions. At 
each such transition the winding number susceptibility is 
enhanced. 


ground-state phase diagram was discussed in Ref. [29|. 

As we have already mentioned in the discussion of the 
phase diagram, the system goes through different mag¬ 
netic configurations upon lowering the field through the 
skyrmion phase to zero. The precise configuration de¬ 
pends on the ratio of system size to pitch length. In 
Fig. [3] we show the situation for a wire geometry with 
linear size L = 16 and one periodic boundary at zero 
temperature and for zero anisotropy. Starting from large 
fields (and hence from the polarized phase), at a cer- 
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Figure 3. Magnetization (top row), energy density pE (middle row), and topological charge density pw (bottom row) for 
different values of the field for a square system with linear size L = 16 with periodic boundaries in the horizontal direction. 
The pitch lengh is p = 7 and the temperature is zero. The color coding indicates out of plane magnetization (top row) and 
energy and topological charge densities (middle and bottom row). 


tain field strength (see Fig. [T]) the magnetic configu¬ 
ration changes from being uniformly polarized, into a 
configuration with skyrmions, upon lowering the field. 
As the distance between skyrmions depends on the field 
(with larger field leading to larger distance) the num¬ 
ber of skyrmions increases discontinuously as the field 
is lowered. Roughly speaking, the number of skyrmions 
jumps once the field (and thus preferred skyrmion dis¬ 
tance) is low enough to accommodate more skyrmions in 
the confined geometry. Note that this clearly depends 
on system size. For lower fields, the skyrmion configu¬ 
ration changes into a spiral state via a state where the 
skyrmions are elongated. Such extended skyrmions are 
reminiscent of “fingers” in liquid crystals that appear in 
many varieties^^— Also note that for low enough fields 
half skyrmions appear at the edge of the system. For 
very low fields, the spiral state is stabilized and in the 
middle of the system the elongated skyrmions and spi¬ 
rals orient themselves 90 degrees with respect to the lat¬ 
tice to maximize the period of the spiral (and thus min¬ 
imize exchange energy). Similar anisotropies may very 
well be present in some materials and in our simulations 
result from the underlying lattice. For simulations that 
would need to give quantitative predictions for contin¬ 
uum systems, one could add additional terms that make 
the exchange interactions more isotropicHere, we do 
not pursue this route as we are interested in the qualita¬ 
tive features of the phases and cascade of phase transi¬ 


tions. At the edge of the wire the influence of exchange is 
less important with respect to DM interactions and the 
spirals are parallel to the edge. Because of the strong 
easy-plane anisotropy of MnSi thin films, Wilson et al. 
considered helicoids with wave vector perpendicular to 
the edges, and, moreover, only considered edge states 
without internal structureIn the geometry we con¬ 
sider, the edge states may also consist of half skyrmions, 
and, as a result, the spiral states that finger out of the 
skyrmions and half skyrmions upon lowering the field will 
have their wave vector parallel to the edges of the wire. 
Again, we note that several of the features we dicuss, such 
as half and elongated skyrmions, and spiral orientation 
at the edges, were very recently experimentally observed 
by Du et al. for Bloch skyrmions^^ 

In Fig. [3] we show, in addition to the magnetization in 
the top row, the energy density pe [the expectation value 
of the summand in Eq. (|4])] in the middle row of figures, 
and the topological charge (winding number) density p^ 
[integrand in Eq. ([5])] in the bottom row, where the 
brighter colors indicate higher values. At the skyrmion 
(and half skyrmion) positions the energy is largest as the 
spin at the skyrmion core points opposite to the exter¬ 
nal field. We also note that the figure of the magne¬ 
tization clearly shows edge states where the magnetiza¬ 
tion tilts away from the field direction at the boundaries 
of the system^^ The middle row shows that the energy 
density is minimized at the edges. Note that for low 
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enough field the half skyrmions at the edge are formed 
from these edge states, as was also discussed in Ref. |4l| 
for Bloch skyrmions. We note that in structures where 
half skyrmions appear the edges show alternating positive 
and negative contributions to the winding number, while 
in the skyrmion phase the contributions to the winding 
number are positive and due to the skyrmions only. 


IV. DEGENERACIES AND FLUCTUATIONS 

At the transitions between different magnetic configu¬ 
rations as a function of field, two magnetic configurations 
are degenerate. Such degeneracies occur generically in a 
confining geometry as a function of field, ratio of DM 
to exchange interaction, or system size, because the ge¬ 
ometry prevents the skyrmions from reaching their pre¬ 
ferred (in an infinitely large system) position. This leads, 
at some particular set of values of the parameters, e.g. 
to a degeneracy between a state with fewer but larger 
skyrmions and a state with more but smaller skyrmions. 
An example of two degenerate spin configurations is given 
in Fig. m In this figure, the skyrmions are elliptical 
as the edge states effectively push the skyrmions away 
from the boundaries^ thereby deforming the skyrmions. 
Skyrmion deformation due to anisotropy was discussed 
in Ref. HU We now turn to the question of how fluc¬ 
tuations lift zero-temperature degeneracies between two 
magnetic configurations, in particular for degenerate con¬ 
figurations containing a different number of skyrmions. 

To investigate this in detail we look again at a wire 
geometry of size L x L = 16 x 16 with KJjD^ = 0.0 and 
BJjD^ = 0.5. We vary the pitch length p and determine 
the classical energy of different skyrmion configurations. 
These energies are found by initiating the simulations 
with a certain number of skyrmions and using simulated 
annealing to force the configuration to an energetic lo¬ 
cal minimum at zero temperature. The upper panel of 
Fig. [U (a) shows the classical ground-state energies for 
various configurations as a function of the pitch length 
relative to the energy of a system without any skyrmions. 
There is a range of pitch lengths 5.5 13.5 for which 

the classical ground state is a configuration containing 
skyrmions. 


We note that a similar result can be obtained by vary¬ 
ing the system size rather than the pitch length, since the 
ratio between these two length scales determines the pre¬ 
ferred amount of skyrmions in the system. This means 
that for a given material with certain values for J, D, and 
K and some applied magnetic field such that skyrmions 
are expected, wires can be made with a certain thickness 
so that their classical magnetic ground state is degener¬ 
ate. 

To investigate the effect of quantum and thermal spin- 
wave (magnon) fluctuations at zero and nonzero tem¬ 
perature, respectively, we use the method outlined in 




* 'rf” f \ 
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Figure 4. (a): Difference in energy (upper panel: classical 
energy, lower panel: energy including quantum correction) 
between configurations with skyrmions and the state without 
skyrmions at zero temperature versus pitch length, for vari¬ 
ous numbers of skyrmions and system size LxL = 16xl6 
with one periodic and one open boundary and parameters 
KJjD^ = 0.0 and BJ/D^ = 0.5. The blue circular, yel¬ 
low rectangular, green diamond-shaped, and red triangu¬ 
lar data points correspond to systems with 0, 1, 2, and 4 
skyrmions respectively, (b) and (c): Magnetic configura¬ 
tion of two classically-degenerate ground states with energy 
E = —596.17J of a system at pitch length p = 9.336 deep in 
the skyrmion regime with one periodic (horizontal direction) 
and one open (vertical direction) boundary. 


Refs. 0 and 0. We quantize the spins and use a 
Holstein-Primakoff transformation to bosonic operators 
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cly and a].. It is given by 

Sp * - S ITLy 5 

S~ = dr\/25' — Uy , (6) 

where S~ is the usual spin-lowering operator and hp = 
djdp. Moreover, rip denotes the classical spin configu¬ 
ration that is found from the simulations at zero tem¬ 
perature and S is the spin quantum number (which we 
take equal to one as the simulations are done for nor¬ 
malized spins). We insert the above transformation in 
the hamiltonian and keep terms up to quadratic order in 
the creation and annihilations operators which amounts 
to a linear approximation in which interactions between 
spin waves are neglected, which is sufficient for low tem¬ 
peratures. The hamiltonian acquires terms ^ aa and 
^ a) a) that are removed by a Bogoliubov transforma¬ 
tion to new bosonic operators yj and that respectively 
create and annihilate a spin wave with energy e^. Here, 
i is an index that labels the spin-wave modes. After 
the Bogoliubov transformation, the hamiltonian is in the 
above-mentioned harmonic approximation given hy^ 

if = i;ei + ^o + Ee»(%T + I) • (7) 

i ^ '' 

The first term is the classical ground-state energy as 
found in the simulations whereas the second term is a 
quantum contribution that arises in the Bogoliubov ap¬ 
proach outlined aboved. This latter contribution is ab¬ 
sent at the classical level. 

Using the above hamiltonian, the ground-state energy 
including quantum corrections is found to be 

E = E,i + Eo + Y,^ ■ ( 8 ) 

i 

At a point where two magnetic configurations in our sim¬ 
ulations are found to be degenerate, the first term Ec\ is 
equal for them. The spin-wave spectrum (and hence the 
quantum correction represented by the last two terms in 
the above expression for the energy) is, however, gener- 
ically different for two classically degenerate configura¬ 
tions, so that quantum fluctuations may indeed remove 
classical degeneracies. In the lower panel of Fig. [4| (a) 
we show the energy including quantum corrections for 
magnetic configurations containing up to four skyrmions 
as a function of pitch length and for the same parame¬ 
ters as the upper panel of this figure. This result shows 
that the energies are shifted by the quantum corrections. 
Moreover, they are shifted in such a way that the region 
where the configuration with four skyrmions is the true 
lowest-energy state is enlarged with respect to the clas¬ 
sical result. This conclusion is in line with the finding 


At nonzero temperatures, the probability of the sys- 


of Roldan-Molina et al. (Hi that quantum fluctuations 
stabilize skyrmion textures over the collinear ferromag¬ 
netic phase. Loosely speaking this comes about because 
larger gradients in the spin direction lead to more quan¬ 
tum fluctuations. 

Having discussed the effect of quantum, i.e., zero- 
temperature, fluctuations, we now turn to thermal fluc¬ 
tuations. To investigate which configuration is entrop- 
ically preferred, we i) compute the entropy due to the 
spin waves around two degenerate skyrmion configura¬ 
tions (an approach appropriate at low temperatures), and 
ii) directly measure the relative probability of occurrence 
of configurations at nonzero temperatures within our sim¬ 
ulations (and approach valid at intermediate and high 
temperatures). In the first approach the spin waves (or 
magnons) are considered as non-interacting bosonic par¬ 
ticles thus taking into account their quantum statistics. 
We refer to this approach as the spin-wave analysis. 

We look at a system with p = 9.336 and p = 7.284. At 
these values of the pitch length, the classically degenerate 
ground states contain either one or two skyrmions for p = 
9.336 as in Fig. [4| or two or four skyrmions for p = 7.284. 
At nonzero temperature, the total energy of the system 
is in the harmonic approximation given by 

E(T) = -Eel + ^i'^i 5 (9) 

with the Bose-Einstein distribution function rii = 
[exp(/3ei) — 1]“^. In the above equation and following 
discussion we neglect the quantum corrections to the en¬ 
ergy discussed previously. The spin-wave entropy is given 
by 

^(r) =-fcs y] [ni In ni - (1 + ni)ln(l + ni)] . (10) 

i 

Since the simulations are classical we expect them to 
correspond to the Rayleigh-Jeans limit, Ui of 

the above formulas, which leads to equipartition of en¬ 
ergy such that in this limit the total energy is equal to 
NksT, with N the number of spins. In Fig.[5]we display 
the average energies of the degenerate ground states as 
a function of temperature. As expected, at low temper¬ 
atures the simulation agrees with this classical equipar¬ 
tition result, whereas the quantum-mechanical result is 
suppressed with respect to equipartition. This is because 
apart from a few (nearly) zero modes, most spin-wave ex¬ 
citations have energies ^ E, such that at temperatures 
ksT <C B our simulations overestimate their contribu¬ 
tion to the energy. The deviation of the simulations from 
equipartition at high temperatures is because in this limit 
the harmonic approximation starts to break down. 


tern being in one out of two degenerate ground states 
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Figure 5. The energyof various spin configurations, with pitch length p = 9.336 in (a) and p = 7.284 in (b), as a function of 
temperature ksT/ J and where Ed is put equal to zero. In both cases the parameters are KJjD^ = 0.0 and BJjD^ = 0.5. The 
data points represent data from Monte Carlo simulations that are in agreement with results from the equipartition theorem at 
low temperatures (shown by the black dashed line). The solid black lines represent the energies resulting from the spin-wave 
analysis. 



Figure 6. Ratio of probabilities of having a certain amount of skyrmions for two different situations with parameters KJjD^ — 
0.0 and BJ/D^ = 0.5 and pitch length p = 9.336 for (a) and p = 7.284 for (b) as a function of temperature. The dots correspond 
to results from the Monte Carlo simulations, the solid-line from the spin-wave analysis not including translation entropy, and 
the dashed line to the zero-temperature classical limit. 


depends on their difference in entropy. Besides the en¬ 
tropy due to magnons there is also translational entropy 
that depends on the number of skyrmions, and that needs 
to be included in the overall entropy of a configuration. 
The simulations automatically include this extra entropy, 
but it is not accounted for in the expression in Eq. m 
and needs to be included on top of this expression. For 
example, the skyrmion configurations in Fig. [4| can be 
translated in the periodic direction by a lattice constant. 
For a single skyrmion there are 16 (for system size 16 
in the periodic direction) translations, whereas for two 
skyrmions the configuration can be mapped onto itself 
by a translation over 8 lattice site in combination with 


a reflection. Hence, the configurations with one and two 
skyrmions have the same translational entropy. Similar 
counting leads to the conclusion that there are twice as 
much translations possible for the situation of two, as 
compared to four skyrmions so that the translational en¬ 
tropy of the configuration with two skyrmions is ln2 
higher than that of four skyrmions. 

Fig. [6] displays the ratio of probabilities for occurence 
of skyrmion configurations with different number of 
skyrmions: two and one skyrmions for p = 9.336, and 
four and two skyrmions for p = 7.284, as a function 
of temperature. In this figure, both simulation results 
(dots) and results from the spin-wave analysis without 
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the translational entropy are shown (black solid line). 
The probability that results from the classical limit of 
the entropy [Eq. (p!Ul) with the replacement rii 
dashed black line] leads to a probability ratio that agrees 
with the low-temperature limit of our simulations, where 
the ratio was measured by cooling a system down 10^ 
times while measuring the number of skyrmions in the 
system. According to the spin-wave analysis, the proba¬ 
bility ratio is in the low-temperature only determined by 
translation entropy. Both the quantum-mechanical result 
for the entropy as well as the simulations show a peak 
around the ferromagnetic phase transition ksT/ J ^ 1. 
At moderate temperatures (O.IJ —0.5J) our results show 
that the configuration with the most skyrmions is entrop- 
ically favored. 



Figure 7. Winding number (left axis) and average magnetiza¬ 
tion in the ^-direction (right axis) as a function of magnetic 
field for the same parameters as Fig. O Labels a-f refer to 
states in Fig.[3l All parameters are taken the same as for the 
results Fig. [3l 


V. CONCLUSIONS, DISCUSSION, AND 
OUTLOOK 

In this article we have shown that in finite systems 
and at elevated temperatures skyrmions are present in a 
large part of the phase diagram. We have also discussed 
how the magnetic field tunes the system through a cas¬ 
cade of transitions between different magnetic configu¬ 
rations, and how zero-temperature degeneracies between 
such magnetic configurations are lifted by fluctuations. 
Throughout this article we have focused on PMA ma¬ 
terials where the DM interactions are believed to arise 
due to interfaces between very thin layers of magnetic 
materials (such as Co) and materials with strong spin- 
orbit couping (such as Pt)^i2i“— Such DM interactions 
give rise to Neel skyrmions. Because the magnetic lay¬ 
ers in these system are very thin (only a few atoms), our 
two-dimensional treatment is appropriate. In particular. 


the two-dimensional nature and form of the DM inter¬ 
actions prevent the formation of conical phases in the 
phase diagram. This situation is different from the situ¬ 
ation of thin films of MnSi in which Bloch skyrmions are 
stabilized and conical phases are present 

In this article we have focused on the confined geom¬ 
etry of a wire. The cascade of transitions as a function 
of field is generic and appears for any confined geome¬ 
try because the field influences the preferred skyrmion 
distance. Which particular magnetic structures occur in 
the cascade depends on the confined geometry, however, 
and on the form of the DM interactions and anisotropy. 
In particular, we expect that the half-skyrmions that we 
have found will only appear in a wire geometry. 

Jumps in the magnetoresistance in MnSi nanowires 
have been observed recently, and were attributed to 
changes in the number of skyrmions in the magnetic 
configurationMeasurements of changes in the topolog¬ 
ical Hall effect due to the appearance of single skyrmions 
in a confined geometry were performed on FeGe4i Moti¬ 
vated by these experiments we now investigate whether 
the transitions between various configurations that ap¬ 
pear in our system as a function of field (see Fig. n can 
be detected electrically. To this end, we compute the to¬ 
tal winding number and the total magnetization in the 
2 ;-direction as a function of field for the same parame¬ 
ters as Fig. [3] (see Fig. [7j). The total winding number 
determines (up to prefactors) the topological Hall sig¬ 
nal, provided spin-orbit coupling is smal L^^d^ The result 
in Fig. [3 clearly shows jumps in the winding number as 
the magnetic configuration undergoes a structural tran¬ 
sition. For PMA materials, however, the topological Hall 
signal is expected to be very small4i An alternative for 
electrical detection of the magnetization configuration is 
then the anomalous Hall signal that is proportional to 
the total magnetization in the ^-direction for our geom¬ 
etry. Fig. [7| shows that this quantity also jumps as the 
magnetization configuration undergoes a transition as a 
function of field. Based on this analysis we conclude that 
the transitions between various magnetic configurations 
that have we found may be detected electrically. 

Using experimental parameters for PMA materials 
from Ref. [13 we estimate that ^ —1 and 

BJjD^ ^ 10“^ — 10“^ for these experiments. The 
route to observe skyrmions in these systems would there¬ 
fore be to increase the field and lower the anisotropy 
(or preferrably make it easy-plane). Very recently, 
Moreau-Luchaire et al. have reported the observation 
of skyrmions at room temperature in multilayers of Co 
and Pt with PMA4^ The skyrmions observed in these 
measurements are rather large and stabilized as a result 
of both dipole-dipole and DM interactions, and are there¬ 
fore in a somewhat different regime from the skyrmions 
that we have studied in this article. 

In future work, we will investigate how current-induced 
torques manipulate the skyrmionic magnetic structures 
we have found. Finally, motivated by the recent exper¬ 
imental results of Du et al^ we also intend to consider 
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Bloch skyrmions. 
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